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Abstract 

In this paper we review recent developments in the statistical theory of weakly nonlinear dispersive waves, the subject 
known as Wave Turbulence (WT). We revise WT theory using a generalisation of the random phase approximation 
(RPA). This generalisation takes into account that not only the phases but also the amplitudes of the wave Fourier 
modes are random quantities and it is called the "Random Phase and Amplitude" approach. This approach allows to 
systematically derive the kinetic equation for the energy spectrum from the the Peierls-Brout-Prigogine (PBP) equation 
for the multi-mode probability density function (PDF). The PBP equation was originally derived for the three- wave 
systems and in the present paper we derive a similar equation for the four-wave case. Equation for the multi-mode 
PDF will be used to validate the statistical assumptions about the phase and the amplitude randomness used for 
WT closures. Further, the multi-mode PDF contains a detailed statistical information, beyond spectra, and it finally 
allows to study non-Gaussianity and intermittency in WT, as it will be described in the present paper. In particular, 
we will show that intermittency of stochastic nonlinear waves is related to a flux of probability in the space of wave 
amplitudes. 

1 Introduction 

Imagine surface waves on sea produced by wind of moderate strength, so that the surface is smooth and there is 
no whitecaps. Typically, these waves exhibit great deal of randomness and the theory which aims to describe their 
statistical properties is called Wave Turbulence (WT). More broadly, WT deals the flelds of dispersive waves which 
are engaged in stochastic weakly nonlinear interactions over a wide range of scales in various physical media. Plentiful 
examples of WT are found in oceans, atmospheres, plasmas and Bose- Einstein condensates ITll^l^llllSl 1^171 1101 1111 IT^ 
I30| . WT theory has a long and exciting history which started in 1929 from the pioneering paper of Peierls who derived 
a kinetic equation for phonons in solids |19|. In the 1960's these ideas have been vigorously developed in oceanography 
013111311301 and in plasma physics EHIIlEl- First of aU, both the ocean and the plasmas can support great many 
types of dispersive propagating waves, and these waves play key role in turbulent transport phenomena, particularly 
the wind-wave friction in oceans and the anomalous diffusion and thermo-conductivity in tokamaks. Thus, WT kinetic 
equations where developed and analysed for different types of such waves. A great development in the general WT 
theory was done by in the papers of Zakharov and Filonenko (S). Before this work it was generally understood that the 
nonlinear dispersive waveflelds are statistical, but it was also thought that such a "gas" of stochastic waves is close to 
thermodynamic equilibrium. Zakharov and Filonenko |S] were the flrst to argue that the stochastic wavefields are more 
like Kolmogorov turbulence which is determined by the rate at which energy cascades through scales rather than by a 
thermodynamic "temperature" describing the energy equipartition in the scale space. This picture was substantiated 
by a remarkable discovery of an exact solution to the wave-kinetic equation which describes such Kolmogorov energy 



cascade. These solutions are now commonly known as Kolmogorov-Zakharov (KZ) spectra and they form the nucleus 
of the WT theory. 

Discovery of the KZ spectra was so powerful that it dominated the WT theory for decades thereafter. Such spectra 
were found for a large variety of physical situations, from quantum to astrophysical applications, and a great effort 
was put in their numerical and experimental verification. For a long time, studies of spectra dominated WT literature. 
A detailed account of these works was given in which is the only book so far written on this subject. Work on 
these lines has continued till now and KZ spectra were found in new applications, particularly in astrophysics )lti| . 
ocean interior |18| and even cosmology I27| . However, the spectra do not tell the whole story about the turbulence 
statistics. In particular, they do not tell us if the wavefield statistics is Gaussian or not and, if not, in what way. This 
question is of general importance in the field of Turbulence because it is related with the intermittency phenomenon, 
- an anomalously high probability of large fluctuations. Such "bursts" of turbulent wavefields were predicted based 
on a scaling analysis in 32 and they were linked to formation of coherent structures, such as whitecaps on sea |28| 
or collapses in optical turbulence [7]. To study these problems qualitatively, the kinetic equation description is not 
sufficient and one has to deal directly with the probability density functions (PDF). 

In fact, such a description in terms of the PDF appeared already in the the Peierls 1929 paper simultaneously with 
the kinetic equation for waves |19l . This result was largely forgotten by the WT community because fine statistical 
details and intermittency had not interested turbulence researchers until relatively recently and also because, perhaps, 
this result got in the shade of the KZ spectrum discovery. However, this line of investigation was continued by Brout 
and Prigogine I2U| who derived an evolution equation for the multi-mode PDF commonly known as the Brout-Prigogine 
equation. This approach was applied to the study of randomness underlying the WT closures by Zaslavski and Sagdeev 
|21| . All of these authors, Peierls, Brout and Prigogine and Zaslavski and Sagdeev restricted their consideration to 
the nonlinear interaction arising from the potential energy only (i.e. the interaction Hamiltonian involves coordinates 
but not momenta). This restriction leaves out the capillary water waves, Alfven, internal and Rossby waves, as well 
as many other interesting WT systems. Recently, this restriction was removed by considering the most general three- 
wave Hamiltonian systems |24| . It was shown that the multi-mode PDF still obeys the Peierls-Brout-Prigogine (PBP) 
equation in this general case. This work will be described in the present review. We will also present, for the first time, 
a derivation of the evolution equation for the multi-mode PDF for the general case of four wave-systems. This equation 
is applicable, for example, to WT of the deep water surface gravity waves and waves in Bose-Einstein condensates or 
optical media described by the nonlinear Schroedinger (NLS) equation. 

We will also describe the analysis of papers |24| of the randomness assumptions underlying the statistical WT 
closures. Previous analyses in this field examined validity of the random phase assumption |2UI 121) without devoting 
much attention to the amplitude statistics. Such "asymmetry" arised from a common mis-conception that the phases 
evolve much faster than amplitudes in the system of nonlinear dispersive waves and, therefore, the averaging may be 
made over the phases only "forgetting" that the amplitudes are statistical quantities too (see e.g. 0). This statement 
become less obvious if one takes into account that we are talking not about the linear phases ujt but about the phases 
of the Fourier modes in the interaction representation. Thus, it has to be the nonlinear frequency correction that 
helps randomising the phases |21|. On the other hand, for three- wave systems the period associated with the nonlinear 
frequency correction is of the same order in small nonlinearity e as the nonlinear evolution time and, therefore, phase 
randomisation cannot occur faster that the nonlinear evolution of the amplitudes. One could hope that the situation is 
better for 4-wave systems because the nonlinear frequency correction is still ~ but the nonlinear evolution appears 
only in the order. However, in order to make the asymptotic analysis consistent, such correction has to be removed 
from the interaction-representation amplitudes and the remaining phase and amplitude evolutions are, again, at the 
same time scale (now 1/e*). This picture is confirmed by the numerical simulations of the 4-wave systems |26l I23| 
which indicate that the nonlinear phase evolves at the same timescale as the amplitude. Thus, to proceed theoretically 
one has to start with phases which are already random (or almost random) and hope that this randomness is preserved 
over the nonlinear evolution time. In most of the previous literature prior to )24l such preservation was assumed but 
not proven. Below, we will describe the analysis of the extent to which such an assumption is valid made in 1241 . 

We will also describe the results of '^'2T- who derived the time evolution equation for higher-order moments of the 
Fourier amplitude, and its application to description of statistical wavefields with long correlations and associated 
"noisiness" of the energy spectra characteristic to typical laboratory and numerical experiments. We will also describe 
the results of |23l about the time evolution of the one-mode PDF and their consequences for the intermittency of 
stochastic nonlinear wavefields. In particular, we will discuss the relation between intermittency and the probability 
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fluxes in the amplitude space. 



2 Setting the stage I: Dynamical Equations of motion 

Wave turbulence formulation deals with a many-wave system with dispersion and weak nonlinearity. For systematic 
derivations one needs to start from Hamiltonian equation of motion. Here we consider a system of weakly interacting 
waves in a periodic box 

where c; is often called the field variable. It represents the amplitude of the interacting plane wave. The Hamiltonian 
is represented as an expansion in powers of small amplitude, 

H = H2 + Ha+H4+H5 + (2) 
where Hj is a term proportional to product of j amplitudes ci, 

oo 

n ,12, 13 ^ ■■ -In , PI ,P2 ■Pm = i 

where qi, q2, qa, ■ ■ ■ Qn and pi,p2,. . -Pm are wavevectors on a d- dimensional Fourier space lattice. Such general j-wave 
Hamiltonian describe the wave- wave interactions where n waves collide to create m waves. Here T^i'p2 'pl^ represents 
the amplitude of the n ^ m process. In this paper we are going to consider expansions of Hamiltonians up to forth 
order in wave amplitude. 

Under rather general conditions the quadratic part of a Hamiltonian, which correspond to a linear equation of 
motion, can be diagonalised to the form 

oo 

"^2 = ^ tJ„|c„|^. (3) 
n=l 

This form of Hamiltonian correspond to noninteracting (linear) waves. First correction to the quadratic Hamiltonian 
is a cubic Hamiltonian, which describes the processes of decaying of single wave into two waves or confluence of two 
waves into a single one. Such a Hamiltonian has the form 

oo 

(4) 

where e ^ 1 is a formal parameter corresponding to small nonlinearity ( e is proportional to the small amplitude whereas 
c„ is normalised so that c„ ~ 1.) Most general form of three-wave Hamiltonian would also have terms describing the 
confluence of three waves or spontaneous appearance of three waves out of vacuum. Such a terms would have a form 



Z ,m,n — 1 



It can be shown however that for systems that are dominated by three-wave resonances such terms do not contribute 
to long term dynamics of systems. We therefore choose to omit those terms. 

The most general four-wave Hamiltonian will have 3, 3^1, 2^2, 4^0 and 0^4 terms. Nevertheless 
1^3, 3^1, 4^0 and 0^4 terms can be excluded from Hamiltonian by appropriate canonical transformations, so 
that we limit our consideration to only 2^2 terms of 7^4, namely 

oo 

m,ri,/.i,^'— 1 

(5) 
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It turns out that generically most of the weakly nonlinear systems can be separated into two major classes: the ones 
dominated by three-wave interactions, so that Tiz describes all the relevant dynamics and Tii can be neglected, and 
the systems where the three-wave resonance conditions cannot be satisfied, so that the Tiz can be eliminated from 
a Ifamiltonian by an appropriate near-identical canonical transformation |25|. Consequently, for the purpose of this 
paper we are going to neglect either 7^3 or TLi, and study the case of resonant three- wave or four- wave interactions. 

Examples of three- wave system include the water surface capillary waves, internal waves in the ocean and Rossby 
waves. The most common examples of the four-wave systems are the surface gravity waves and waves in the NLS model 
of nonlinear optical systems and Bose-Einstein condensates. For reference we will give expressions for the frequencies 
and the interaction coefficients corresponding to these examples. 

For the capillary waves we have 



and 



where 



V' = 

' m.n 



8nV2a 



Nl/2 



Lkr, 



(^kj^kn^^^^kl (klkra)^^^kn {klkn^^^^k-n 



— (km ■ kri) kjnkn 



and a is the surface tension coefficient. 
For the Rossby waves |13II14| . 



and 



-^\k, k 

47r' 



knx I ^ X 



mx r\'nx 



1 +p2fc2 



where /3 is the gradient of the Coriolis parameter and p is the Rossby deformation radius. 



The simplest expressions correspond to the NLS waves |15l 171. 



VF'" = 1. 



The surface gravity waves are on the other extreme. The frequency is = 
notoriously long expressions which can be found in [Tl I17|. 

2.1 Three- wave case 

When Ti — Ti.2 + "Hz we have Hamiltonian in a form 



(6) 
(7) 
(8) 

(9) 
(10) 

(11) 

gk but the matrix element is given by 



H = ^a;„|c„|^ + e ^ \V^nClCmCn5m+n + C.C. 
n=\ I ,m,,n—l 

Equation of motion ici = is mostly conveniently represented in the interaction representation, 

oc 



(12) 



where aj = CjC^'^^* is the complex wave amplitude in the interaction representation, l,m,n G Z are the indices 
numbering the wavevectors, e.g. km = 2-Km/L, L is the box side length, ujln„ = — ^^km ~ '^fem '^i ~ is 
the wave linear dispersion relation. Here, V^^i ~ 1 is an interaction coefficient and e is introduced as a formal small 
nonlinearity parameter. 
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2.2 Four- wave case 

Consider a weakly nonlinear wavefield dominated by the 4-wave interactions, e.g. the water-surface gravity waves 
IHl l29ll^ . Langmuir waves in plasmas Tj'SI and the waves described by the nonlinear Schroedinger equation [7|. The 
a Hamiltonian is given by (in the appropriately chosen variables) as 



n—l m,,rL,(i,i' — l 



(13) 



As in three-wave case the most convenient form of equation of motion is obtained in interaction representation, c; — 
bie-^^'*, so that 

ibi = eJ2 WlZbc^b^b^e^^^'SlZ (14) 



where Wj^Z ~ 1 is an interaction coefficient, tjj," — loi + coa ~ lo^ — uj,y. We are going expand in e and consider the 
long-time behaviour of a wave field, but it will turn out that to do the perturbative expansion in a self-consistent 
manner we have to renormalise the frequency of l|14|l as 



— typ°aaaMai^e'"^""'^^" — Qiai, (15) 



where ai = foje'"'* , = lu''^, + ^1 +Qa - - ^1^, and 

=2e^M//;4 (16) 

is the nonlinear frequency shift arising from self-interactions. 

3 Setting the stage II: Statistical setup 

In this section we are going to introduce statistical objects that shall be used for the description of the wave systems, 
PDF's and a generating functional. 

3.1 Probability Distribution Function. 

Let us consider a wavefield a(x, t) in a periodic cube of with side L and let the Fourier transform of this field be 
ai{t) where index marks the mode with wavenumber fc; = 2-kI/L on the grid in the d- dimensional Fourier space. 

For simplicity let us assume that there is a maximum wavenumber kmax (fixed e.g. by dissipation) so that no modes 
with wavenumbers greater than this maximum value can be excited. In this case, the total number of modes is 
N = {kmax/i^LY . Correspondingly, index I will only take values in a finite box, I £ Bn Z"^ which is centred at 
and all sides of which are equal to kmax/i^L — N^^^ . To consider homogeneous turbulence, the large box limit N ^ 00 
will have to be taken. ^ 

Let us write the complex a; as a; = Aitpi where Ai is a real positive amplitude and tpi is a phase factor which 
takes values on S^, a unit circle centred at zero in the complex plane. Let us define the A^'-mode joint PDF p'^' 
as the probability for the wave intensities Af to be in the range (s;,s; -I- dsi) and for the phase factors i/i; to be on 
the unit-circle segment between ^; and (,1 + d^; for all I £ Bn- In terms of this PDF, taking the averages will involve 
integration over all the real positive s;'s and along all the complex unit circles of all ^;'s. 




^It is easily to extend the analysis to the infinite Fourier space, kmax = 00. In this case, the full joint PDF would still have to be defined 
as a A'^ ^ 00 limit of an TV-mode PDF, but this limit would have to be taken in such a way that both kmax and the density of the Fourier 
modes tend to infinity simultaneously. 
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where notation /{^^, ?/>} means that / depends on aU Af's and all i/'i's in the set {Af,tlji; 1€Bn} (similarly, {s, ^} means 
{si,^i;l £ Bn}, etc). The full PDF that contains the complete statistical information about the wavefield a(x, t) in 
the infinite x-space can be understood as a large-box limit 

lim 

Ai — >oo 

i.e. it is a functional acting on the continuous functions of the wavenumber, Sk and ^fe. In the the large box limit there 
is a path-integral version of 11711 , 

{f{A\^}) = J vsjm (18) 

The full PDF defined above involves all A'^ modes (for either finite A'^ or in the N —* oo limit). By integrating out all 
the arguments except for chosen few, one can have reduced statistical distributions. For example, by integrating over 
all the angles and over all but M amplitudes,we have an "M-mode" amplitude PDF, 

= ( n I ds^w M^-i ] ^'"^i'*'?}' (19) 

which depends only on the M amplitudes marked by labels j'l, J2, • . . jJm&Bn- 



3.2 Definition of an ideal RPA field 

Following the approach of I23| . we now define a "Random Phase and Amplitude" (RPA) field. We say that the 
field a is of RPA type if it possesses the following statistical properties: 

1. All amplitudes Ai and their phase factors ipi are independent random variables, i.e. their joint PDF is equal to 
the product of the one-mode PDF's corresponding to each individual amplitude and phase, 

2. The phase factors tpi are uniformly distributed on the unit circle in the complex plane, i.e. for any mode I 

P/'^)(CO = 1/2^. 

Note that RPA does not fix any shape of the amplitude PDF's and, therefore, can deal with strongly non-Gaussian 
wavefields. Such study of non-Gaussianity and intermittency of WT was presented in |22II23| and will not be repeated 
here. However, we will study some new objects describing statistics of the phase. 

In |22II23| RPA was assumed to hold over the nonlinear time. In |24| this assumption was examined a posteriori, i.e. 
based on the evolution equation for the multi-point PDF obtained with RPA initial fields. Below we will describe this 
work. We will see that RPA fails to hold in its pure form as formulated above but it survives in the leading order so 
that the WT closure built using the RPA is valid. We will also see that independence of the the phase factors is quite 
straightforward, whereas the amplitude independence is subtle. Namely, M amplitudes are independent only up to a 
0{M/N) correction. Based on this knowledge, and leaving justification for later on in this paper, we thus reformulate 
RPA in a weaker form which holds over the nonlinear time and which involves Af-mode PDF's with M <^ N rather 
than the full A^-mode PDF. 

^We keep the same acronym as in related "Random Phase Approximation" but now interpret it differently because (i) we emphasise the 
amplitude randomness and (ii) now RPA is a defined property of the field to be examined and not an approximation. 
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3.3 Definition of an essentially RPA field 

We will say that the field a is of an "essentially RPA" type if: 

1. The phase factors are statistically independent and uniformly distributed variables up to O(e^) corrections, i.e. 

'P'-'^Hs,^} = 7T5rTF^*'^''"{«} [1 + Oie% (20) 

where 




n L i^^'i 1 ^*^'{^'^}' (21) 
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is the A'^-mode amplitude PDF. 

2. The amplitude variables are almost independent is a sense that for each M <^ N modes the M-mode amplitude 
PDF is equal to the product of the one-mode PDF's up to 0{M/N) and o(e^) corrections, 

- ■ • ■ ^t' [1 + 0{M/N) + Oie')]. (22) 

3.4 Why ^p^s and not 0's? 

Importantly, RPA formulation involves independent phase factors ip — e'^ and not phases (j> themselves. Firstly, the 
phases would not be convenient because the mean value of the phases is evolving with the rate equal to the nonlinear 
frequency correction |24|. Thus one could not say that they are "distributed uniformly from —n to tt". Moreover 
the mean fluctuation of the phase distribution is also growing and they quickly spread beyond their initial 27r-wide 
interval |24) . But perhaps even more important, </>'s build mutual correlations on the nonlinear time whereas tp's remain 
independent. Let us give a simple example illustrating how this property is possible due to the fact that correspondence 
between (p and t/i is not a bijection. Let A be a random integer and let ri and r2 be two independent (of N and of 
each other) random numbers with uniform distribution between ~n and tt. Let 

<^i,2 =27rA + ri,2. 

Then 

(,^1.2) =27r{A), 

and 

(,^1,^2) =47r"{A'). 

Thus, 

{<^i<^2) - {<^i){<^2) = 4^'((A^) - (Nf) > 0, 
which means that variables (pi and 02 are correlated. On the other hand, if we introduce 



i>i,2 = e""^-^ , 

then 

(V'1.2) =0, 

and 

= 0. 

{^1^2) -{V'i){V'2) =0, 

which means that variables ipi and -02 are statistically independent. In this illustrative example it is clear that the 
difference in statistical properties between (j) and 1/) arises from the fact that function ^^{(j}) does not have inverse and, 
consequently, the information about N contained in (j> is lost in tp. 

Summarising, statistics of the phase factors tp is simpler and more convenient to use than (p because most of 
the statistical objects depend only on tp. This does not mean, however, that phases (p are not observable and not 
interesting. Phases (p can be "tracked" in numerical simulations continuously, i.e. without making jumps to — tt when 
the phase value exceeds tt. Such continuous in k function (p{k) can achieve a large range of variation in values due to 
the dependence of the nonlinear rotation frequency with k. This kind of function implies fastly fluctuating ^p{k) which 
is the mechanism behind de-correlation of the phase factors at different wavenumbers. 
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3.5 Wavefields with long spatial correlations. 



Often in WT, studies are restricted to wavefields with fastly decaying spatial correlations 2 . For such fields, the 
statistics of the Fourier modes is close to being Gaussian. Indeed, it the correlation length is much smaller than the 
size of the box, then this box can be divided into many smaller boxes, each larger than the correlation length. The 
Fourier transform over the big box will be equal to the sum of the Fourier transforms over the smaller boxes which 
are statistically independent quantities. Therefore, by the Central Limit Theorem, the big-box Fourier transform 
has a Gaussian distribution. Corrections to Gaussianity are small as the ratio of the correlation volume to the box 
volume. On the other hand, in the RPA defined above the amplitude PDF is not specified and can significantly deviate 
from the Rayleigh distribution (corresponding to Gaussian wavefields). Such fields correspond to long correlations 
of order or greater than the box size. In fact, long correlated fields are quite typical for WT because, due to weak 
nonlinearity, wavepackets can propagate over long distance preserving their identity. Moreover, restricting ourselves to 
short-correlated fields would render our study of the PDF evolution meaningless because the later would be fixed at 
the Gaussian state. Note that long correlations modify the usual Wick's rule for the correlator splitting by adding a 
singular cumulant, e.g. for the forth-order correlator, 

where Qc = {At - 2Al). These issues were discussed in detail in 1221. 



3.6 Generating functional. 

Introduction of generating functionals simplifies statistical derivations. It can be defined in several different ways to 
suit a particular technique. For our problem, the most useful form of the generating functional is 

Z(^'{A,M} = y^{n e^'^^^n, (23) 

where {A,/i} = {Xi,^i;l £ Bn} is a set of parameters, Xi£TZ and fiiGZ. 

' {m} leBN 

where {/i} = {fii £ Z\l £ Bn}- This expression can be verified by considering mean of a function using the 

averaging rule and expanding / in the angular harmonics i/;™; m £ Z (basis functions on the unit circle), 

/M^v} = Ef{"^'^} n ^r, (25) 

{m} ieSiv 

where {m} = {m; £ Z;l £ Bn} axe indices enumerating the angular harmonics. Substituting this into II17II with PDF 
given by 1241 and taking into account that any nonzero power of ^; will give zero after the integration over the unit 
circle, one can see that LHS=RHS, i.e. that 12411 is correct. Now we can easily represent 1241 in terms of the generating 
functional, 

= (^'"'KmI n ^r] (26) 

{m} V i68„ / 

where CJ^^ stands for inverse the Laplace transform with respect to all A; parameters and {/i} = {/ii £ Z;l £ Bn} are 
the angular harmonics indices. 

Note that we could have defined Z for all real ^I's in which case obtaining P would involve finding the Mellin 
transform of Z with respect to all /xj's. We will see below however that, given the random-phased initial conditions, Z 
will remain zero for all non-integer ^;'s. More generally, the mean of any quantity which involves a non-integer power 
of a phase factor will also be zero. Expression 12611 can be viewed as a result of the Mellin transform for such a special 
case. It can also be easily checked by considering the mean of a quantity which involves integer powers of tpi's. 
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By definition, in RPA fields all variables Ai and tpi are statistically independent and ipi's are uniformly distributed 
on the unit circle. Such fields imply the following form of the generating functional 

z(^'{A,m} = ^''^'"'{A} n ^(mO, (27) 

where 

Z(^'^^{X} = ( H e'^^') = Z(^'{A,m}|m=o (28) 

is an A^-mode generating function for the amplitude statistics. Here, the Kronecker symbol 5(/i;) ensures independence 
of the PDF from the phase factors tpi. As a first step in validating the RPA property we will have to prove that the 
generating functional remains of form 12711 up to 1/A'' and O(e^) corrections over the nonlinear time provided it has 
this form at i = 0. 



3.7 One-mode statistics 

Of particular interest are one-mode densities which can be conveniently obtained using a one-amplitude generating 
function 

Z(t,A) = {e^l'^'=l'), 

where A is a real parameter. Then PDF of the wave intensities s — \ak\^ at each k can be written as a Laplace 
transform, 

P<--\s,t) = {S{\ak\'-s)) = — Z{X,t)e-'^dX. (29) 

27r Jo 

For the one-point moments of the amplitude we have 

Zx...x\x^o = / s''P^''\s,t)ds, (30) 



where p £ N and subscript A means differentiation with respect to A p times. 

The first of these moments, ~ is the waveaction spectrum. Higher moments measure fiuctuations of 

the waveaction fc-space distributions about their mean values |22|. In particular the r.m.s. value of these fiuctuations 
is 



a, = JMf-n2. (31) 



4 Separation of timescales: general idea 

When the wave amplitudes are small, the nonlinearity is weak and the wave periods, determined by the linear dynamics, 
are much smaller than the characteristic time at which different wave modes exchange energy. In the other words, 
weak nonlinearity results in a timescale separation and our goal will be to describe the slowly changing wave statistics 
by averaging over the fast linear oscillations. To filter out fast oscillations, we will seek seek for the solution at time 
T such that 2% /uj <^ T <^ Tni- Here r„i is the characteristic time of nonlinear evolution which, as we will see later is 
~ l/uje^ for the three- wave systems and ~ l/uje^ for the four- wave systems. Solution at t — T can be sought at series 
in small small nonlinearity parameter e, 

ai(T) =a<°' -heaP' +6'ap\ (32) 

Then we are going to iterate the equation of motion 1121 or 11511 to obtain a'"', a^' and ap' by iterations. 
During this analysis the certain integrals of a type 

t 

fit) = J g{t)eMi^i') 
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will play a crucial role. Following ,2, we introduce 



1) 



(33) 

E{x,y)= I A{x-y)e'y'dt. 
Jo 

We will be interested in a long time asymptotics of the above expressions, so the following properties will be useful: 

lim E{0,x) ^T{tvS{x) +iP{-)), 

T->oo X 



e at — 



and 



and 



lim |A(x)|^ = 2tvT5{x). 

T — »oo 



5 Weak nonlinearity expansion: three-wave case. 

Substituting the expansion 1321 in 11211 we get in the zeroth order 

i.e. the zeroth order term is time independent. This corresponds to the fact that the interaction representation wave 
amplitudes are constant in the linear approximation. For simplicity, we will write a|'''(0) = a;, understanding that a 
quantity is taken at T = if its time argument is not mentioned explicitly. 

Here we have taken into account that aj''''(r) — ai and a^^^(O) = 0. 

The first order is given by 

oo 
m,n— 1 

) , (34) 

Here we have taken into account that af^\T) — ai and a^^'(O) = 0. Perform the second iteration, and integrate over 
time to obtain To calculate the second iterate, write 



• • (2) 



= 1 



m , ri — 1 



(35) 

Substitute 1I34I I into (I35I I and integrate over time to obtain 

oo 

+21/;™ {y^^a^a^,a^S'^+^E[-ujZ^,,-^Z] + ^Vj^^arad^a^Elijj^l,, -ujJ!^]5'^_^^J^ , 

(36) 
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where we used a^^''(0) = 0. Hereafter, we drop super-script (0) in expressions like 1351 for brevity of notations. 



6 Weak nonlinearity expansion: four-wave case 

Substituting 13211 in list we get in the zeroth order af\T) — ai{0), and the first iteration of 11511 gives 



(37) 



Iterating one more time we get 



[fla — 2f2„) / re' i^"^ 

aiJ.v y g 



(38) 



7 Asymptotic expansion of the generating functional. 

Let us first obtain an asymptotic weak-nonlinearity expansion for the generating functional Z{X, /i} exploiting 
separation of the linear and nonlinear time scales. To do this, we have to calculate Z at the intermediate time t ■ 
via substituting into it aj (T) from H32|l For the amplitude and phase "ingredients" in Z we have, 



the 

= T 



gjljl —gjlj -r J -r J I _gj J g a2j), 



and 



where 



==Vf''Ml + ^/3i.+^'/32.), 



Ol2j 



(A, + A^Af + hiaf-af + a^) + ^io!f^^f + (af af )^ 



2A 



Substituting expansions 1391 and 14U1 into the expression for Z, we have 

Z{\, M, T} = X{\, /i, T} + X{A, -M, T} 



with 



4^ 



(0)2 



X{A,^,T} = X{A,M,0} + (27r)2^/ e"''"'"''' [eJi + ^'(Ja + Ja + J4 + Js)]) + 0(e*), 

\||I||<iV / . 



^Hereafter we omit superscript (A'^) in the A^-mode objects if it does not lead to a confusion. 



(39) 

(40) 

(41) 
(42) 
(43) 

(44) 

(45) 
(46) 
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where 
Jl 

J2 

is 



J5 



(0)12'' J' 



2i (0)|2 



_i^^|^{l)|2 



(0)12^' J' 



2|a-l^ 



J 'if} 



^ + 



2 4|a(.'"|4 V 2 



^-1 + 



2|af P 



(a 



(47) 
(48) 
(49) 
(50) 



4 a 



where (•)^ and {•)^ denote the averaging over the initial amplitudes and initial phases (which can be done indepen- 
dently). Note that so far our calculation for Z{T) is the same for the three- wave and for the four- wave cases. Now 
we have to substitute expressions for a'^' and a'^' which are different for the three-wave and the four-wave cases and 
given by linnj and JSEJ respectively. 



8 Evolution of statistics of three-wave systems 
8.1 Equation for the generating functional 

Let us consider the initial fields afc(O) — af^ which are of the RPA type as defined above. We will perform averaging 
over the statistics of the initial fields in order to obtain an evolution equations, first for Z and then for the multi-mode 
PDF. Let us introduce a graphical classification of the above terms which will allow us to simplify the statistical 
averaging and to understand which terms are dominant. We will only consider here contributions from Ji and J2 
which will allow us to understand the basic method. Calculation of the rest of the terms, J3, J4 and J5, follows the 
same principles and can be found in |24| . First, The linear in e terms are represented by Ji which, upon using 1)34^ . 
becomes 

= (n E (A. + ^) (v;i„a™anA^_5^+„ + 2V^^a„a„Af„57+„) a A . (52) 

\ t j,m.n J I ^ 

Let us introduce some graphical notations for a simple classification of different contributions to this and to other 
(more lengthy) formulae that will follow. Combination V^nS-^^^^ will be marked by a vertex joining three lines with 
in-coming j and out-coming m and n directions. Complex conjugate V^„S^^^ will be drawn by the same vertex but 
with the opposite in-coming and out-coming directions. Presence of Uj and aj will be indicated by dashed lines pointing 
away and toward the vertex respectively. ■* Thus, the two terms in formula 1521 can be schematically represented as 
follows. 



*This technique provides a useful classification method but not a complete mathematical description of the terms involved. 



12 




Let us average over all the independent phase factors in the set {'tp}. Such averaging takes into account the statistical 
independence and uniform distribution of variables 'tp. In particular, (ip) = 0, {ipitpm} = and (ipi'ipm) = <5("- Further, 

the products that involve odd number of ?/)'s are always zero, and among the oven products only those can survive 
that have equal numbers of tp's and tp's. These tp's and tp's must cancel each other which is possible if their indices are 
matched in a pairwise way similarly to the Wick's theorem. The difference with the standard Wick, however, is that 
there exists possibility of not only internal (with respect to the sum) matchings but also external ones with tp's in the 
pre-factor Ilipi'' . 

Obviously, non-zero contributions can only arise for terms in which all xp's cancel out either via internal mutual 
couplings within the sum or via their external couplings to the tp's in the Z-product. The internal couplings will indicate 
by joining the dashed lines into loops whereas the external matching will be shown as a dashed line pinned by a blob 
at the end. The number of blobs in a particular graph will be called the valence of this graph. 

Note that there will be no contribution from the internal couplings between the incoming and the out-coming lines 
of the same vertex because, due to the J-symbol, one of the wavenumbers is in this case, which means ^ that V = 0. 
For Ji we have 

Jl = (C'l)v> + (C2)-il>, 

with 



> 



• j + mm 



>__. 



^ •2m 



and 



m 
n 



{C2)v< = 



✓ ■ 



• 3 + < 



• 2n 



which correspond to the following expressions, 

®In the present paper we consider only spatially homogeneous wave turbulence fields, hi spatially homogeneous fields, due to momentum 
conservation, there is no coupling to the zero mode fc = because such coupling would violate momentum conservation. Therefore if one of 
the arguments of the interaction matrix element V is equal to zero, the matrix element is identically zero. That is to say that for any spatially 
homogeneous wave turbulence system Vj^^^^ = VI^_^_q^^ = V^^i^^—q = 0. 



13 



(Ci)v. = Yl {^j + ^)Vln Am An AjAi,J'^+J{fim + l)5{ii„ + 1)5(11,-1) II Silll) 

+ ^(A2m + )VmmAmA2mA'^S{fJ.m + 2)(5(^2m - 1) Yl '^(/^' ) (^3) 

m 2m l5.^m,2m 



and 



+ 2^(A„ + ^)K';'A2„^^A^„5(M2n + l)5(Mn-2) n ^(/^O- (54) 

Because of the 5-symbols involving fi's, it takes very special combinations of the arguments /i in Z{^} for the terms in 
the above expressions to be non-zero. For example, a particular term in the first sum of I53II may be non-zero if two 
/i's in the set {fi} are equal to 1 whereas the rest of them are 0. But in this case there is only one other term in this 
sum (corresponding to the exchange of values of n and j) that may be non-zero too. In fact, only utmost two terms 
in the both 15311 and 15411 can be non-zero simultaneously. In the other words, each external pinning of the dashed 
line removes summation in one index and, since all the indices are pinned in the above diagrams, we are left with no 
summation at all in Ji i.e. the number of terms in Ji is 0(1) with respect to large A'^. We will see later that the 
dominant contributions have 0{N^) terms. Although these terms come in the order, they will be much greater that 
the terms because the limit N ^ oo must always be taken before e ^ 0. 

Let us consider the first of the e^-terms, J2. Substituting 13411 into 14811 . we have 



2 



where 



I j ,7n.rL,K.,iy 

= {Bi + B2 + B2 + B3}^, (55) 



K m K m K 



Bi= >•-*•••; B2= >--^--N and S3 = (56) 



Here the graphical notation for the interaction coefficients V and the amplitude a is the same as introduced in the 
previous section and the dotted line with index j indicates that there is a summation over j but there is no amplitude 
Oj in the corresponding expression. 

Let us now perform the phase averaging which corresponds to the internal and external couplings of the dashed 
lines. For {Bi),p we have 
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m K m K 

/ \ / / '"^ 

• • " ^ m 

nun V 

m 

+ 2 {...^...\ , (57) 



where 



m K 

V i / 1 ^2 



1 2 



We have not written out the third term in H57^ because it is just a complex conjugate of the second one. Observe 
that all the diagrams in the first line of 15711 are 0(1) with respect to large A'^ because all of the summations are lost 
due to the external couplings (compare with the previous section). On the other hand, the diagram in the second line 
contains two purely-internal couplings and is therefore 0{N^). This is because the number of indices over which the 
summation survives is equal to the number of purely internal couplings. Thus, the zero-valent graphs are dominant 
and we can write 

(BiU = n ^(/^O E (^J- + ^j^',)\V;Luf\Ai,ufSi,+„AlAl[l + 0{1/N^)] (58) 

I j,m,,7i 
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Now we are prepared to understand a general rule: the dominant contribution always comes from the graphs with 
minimal valence, because each external pinning reduces the number of summations by one. The minimal valence of 
the graphs in (-62)^ is one and, therefore, {62)^ is order A'^ times smaller than (Bi)^. On the other hand, {B^}^ is of 
the same order as (Bi)^, because it contains zero- valence graphs. We have 

m 



J 



(S3), 



n 

I j,m,n 

Summarising, we have 

J2 =l[^M E + ^?^?) [\Vi„f\^inn\'5U„ + 2H/^nA7„ AlAl [1 + 0{1/N)]. 



(59) 



(60) 



Thus, we considered in detail the different terms involved in J2 and we found that the dominant contributions come 
from the zero-valent graphs because they have more summation indices involved. This turns out to be the general rule 
in both three-wave and four-wave cases and it allows one to simplify calculation by discarding a significant number of 
graphs with non-zero valence. Calculation of terms J3 to J5 can be found in 24^ and here we only present the result, 

Js-^U^M E [-|l^^nl'£(0,a.;^J5^+„A?. + iV;™i^i5(0,~a;7„)57+„(A^--^^)] A^^ x [1 + 0(1/7V)], (61) 



J4 = 0(1) (i.e. J4 is order A'^'^ times smaller than J2 or J3) and 

J5 = 2n<5(M0 E ^,^>^[-mj'K+J^U' + \Vrk\'S7+,\A],\']A',AlAl [l + 0(l/iV)]. 



(62) 



Using our results for Ji — J5 in 14611 and 14511 we have 

Z{T) - Z{0) = E (^^- + ^?^) [\V^r^\'\^Lr^fSUn + 2| ^ A™, j ^<5f+„] 



n 



+4e^ E 



d 



'dX 



( ' 




\dX^ 


dXn) 



dZjO) 
dX, 



a^z(o) 

dXjdXndXk 



+ cc. 



(63) 



Here partial derivatives with respect to A; appeared because of the Ai factors. This expression is valid up to 0(e*) 
and 0{£^ /N) corrections. Note that we still have not used any assumption about the statistics of ^'s. Let us now 
N ^ ca limit followed by T ~ 1/e ^ 00 (we re- iterate that this order of the limits is essential). Taking into account 
that lim E{0,x) = T{tvS{x) + iP{^)), and lim \A{x)\'^ = 2ttTS{x) and, replacing {Z{T) - Z{0))/T by Z we have 



Z = J {{X,+X^~)[\Vl„f5{u^i,„)Si^^^ + 2\Vj::f5(L.Tn)5r^ 



S'Z 



'+"\ SXmSXn 



+2A, 



+ 2X.jXm [-2|V'^„|^(5^ + „(5(t<;^„) + | V^m|^5"+m<5(t^"m) 



s s_ 

SXm 5X n 

6^Z 

SXjSXnSXm 



sz 



(64) 
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Here variational derivatives appeared instead of partial derivatives because of the N ^ oo limit. 

Now we can observe that the evolution equation for Z does not involve fi which means that if initial Z contains factor 
Yli S{fii) it will be preserved at all time. This, in turn, means that the phase factors {tf}} remain a set of statistically 
independent (of each each other and of A's) variables uniformly distributed on S^. This is true with accuracy O(e^) 
(assuming that the A^-limit is taken first, i.e. ^ e^) and this proves persistence of the first of the "essential RPA" 
properties. Similar result for a special class of three-wave systems arising in the solid state physics was previously 
obtained by Brout and Prigogine |2()| . This result is interesting because it has been obtained without any assumptions 
on the statistics of the amplitudes {A} and, therefore, it is valid beyond the RPA approach. It may appear useful in 
future for study of fields with random phases but correlated amplitudes. 



8.2 Equation for the multi-mode PDF 

Taking the inverse Laplace transform of 16411 we have the following equation for the PDF, 
where jFj is a flux of probability in the space of the amplitude sj , 



(65) 



47re2s ■ 



mnl <5(lJ^„„)5: 



5V_ 



y-n 1 2 



j + m 



+2(11/, 



(66) 



This equation is identical to the one originally obtained by Peierls 1191 and later rediscovered by Brout and Prigogine 
|20| in the context of the physics of anharmonic crystals. 



V 



(67) 



(68) 



where 

'S] _ _S 5 5_ 

(5s J g SSj SSm Ss„ 

Zaslavski and Sagdeev |21) were the first to study this equation in the WT context. However, the analysis of |19ll2UlET| 
was restricted to the interaction Hamiltonians of the "potential energy" type, i.e. the ones that involve only the 
coordinates but not the momenta. This restriction leaves aside a great many important WT systems, e.g. the capillary, 
Rossby, internal and MHD waves. Our result above indicates that the Peierls equation is also valid in the most general 
case of 3-wave systems. Note that Peierls form 16711 - 1681 looks somewhat more elegant and symmetric than 1651 - 
1661 . However, form 1651 - 1661 has advantage because it is in a continuity equation form. Particularly for steady state 
solutions, one can immediately integrate it once and obtain. 



where A is an arbitrary functional of s{k) and e;„ 



is the antisymmetric tensor. 



(69) 



1 for ki > km > kn and cyclic permutations of l,Tn,n, 

-1 for all other fc; 7^ km 7^ k„, 
if ki = km^ ki = kn or k^ ~ km- 



In the other words, probability flux F can be an arbitrary solenoidal field in the functional space of s(fe). One can 
see that 1691 is a first order equation with respect to the s-derivative. Special cases of the steady solutions are the 
zero-flux and the constant-flux solutions which, as we will see later correspond to a Gaussian and intermittent wave 
turbulence respectively. 
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Here we should again emphasise importance of the taken order of limits, N ^ <x first and e ^ second. Physically 
this means that the frequency resonance is broad enough to cover great many modes. Some authors, e.g. 19 20 21|. 
leave the sum notation in the PDF equation even after the e ^ limit taken giving (5(tj"„). One has to be careful 
interpreting such formulae because formally the RHS is nill in most of the cases because there may be no exact 
resonances between the discrete k modes (as it is the case, e.g. for the capillary waves). In real finite-size physical 
systems, this condition means that the wave amplitudes, although small, should not be too small so that the frequency 
broadening is sufficient to allow the resonant interactions. Our functional integral notation is meant to indicate that 
the N ^ cx limit has already been taken. 



9 Evolution of statistics of four-wave systems 

In this section we are going to calculate the four-wave analog of PBP equation. These calculations are similar in spirit 
to those presented in the orevious section, but with the four- wave equation of motion p5|l . The calculation will be 
slightly more involved due to the frequency renormalisation, but it will be significantly simplified by our knowledge that 
(for the same reason as for three-wave case) all fi dependent terms will not contribute to the final result. Indeed, any 
term with nonzero /i would have less summations and therefore will be of lower order. On the diagrammatic language 
that would mean that all diagrams with non-zero valence can be discarded, as it was explained in the previous section. 
Below, we will only keep the leading order terms which are ~ A'^^ and wich correspond to the zero-valence diagrams. 
Thus, to calculate Z we start with the Ji — J5 terms 1471 . 148( l . 1491 1 . (ISUl . 15 in in which we put /i = and substutute 
into them the values of a*-^' and o^^' from 13711 and 13811 . For the terms proportional to e we have 

/„{i)^(o)\ _ 

^-2iY, W]^A]Al ■ T + i^j A]t. (70) 

where we have used the fact that A(0) = T. We see that our choice of the frequency renormalisation II6II makes the 
contribution of {af ^a}!^^)^ equal to zero, and this is the main reason for making this choice. Also, this choice of the 
frequency correction simplifies the order and, most importantly, eliminates from them the cumulative growth, 
e.g. 

{.{a'^paff)^ = -2W^iA^T^n, - Ap^jT^ - 4 ^(1V/„"A, A„T)' in the iV 00 limit. 
Finally we get, 

Ji = 0, 

J2 = ^ A,(l + A,.4,^)iiy^^l^lAi:|^5^^A^^^^^„ 
Js = 2J2 X,\W^^\^S';iA^^AliAl-2Al)EiO,Cu';i), 
Ja = 0, 

^5 = Y.^^^^[~'^\^i'^^\^\^k't\'5CAl + \W^^^^^^ (71) 

By putting everything together in 1451 and 1451 1 and explointing the symmeteries introduced by the summation variables, 
we finally obtain in A*' — > cxa limit followed by T ^ 00 

\W^'^f5{Cbii^)5'^^\, — — — — [(A, + A, - A„ - K)Z] dkjdhdk^dK. (72) 

ClAj OAl OAn ClAm 
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By applying to this equation the inverse Laplace transform, we get a four-wave analog of the PBP equation for the 
N-mode PDF: 



■ s ■ 




■ s ■ 


r) 










4 V 




4 / 



V ] dkjdkidkmdkn, 



where 

5_ 5_ 

5sJ^ 5sj 5si Ssm Ss„ 
This equation can be easily written in the continuity equation form I65II with the flux given in this case by 



V dkidkmdkn 



(73) 
(74) 

(75) 



10 Justification of RPA 

Variables Sj do not separate in the above equation for the PDF. Indeed, substituting 

■p(N,a) _ p{a) p(a) p{a) /„„n 



into the discrete version of Ifcifcil we see that it turns into zero on the thermodynamic solution with P^-"^ = ujj exp{—ujjSj). 

However, it is not zero for the one-mode PDF P^-"' corresponding to the cascade-type Kolmogorov-Zakharov (KZ) 

spectrum rij^ , i.e. Pj-"'' = {l/rij") exp{—Sj /rij^) (see the Appendix), nor it is likely to be zero for any other PDF of 
form II75I I. This means that, even initially independent, the amplitudes will correlate with each other at the nonlinear 
time. Does this mean that the existing WT theory, and in particular the kinetic equation, is invalid? 

To answer to this question let us differentiate the discrete version of the equation 1641 with respect to A's to get 
equations for the amplitude moments. We can easily see that 

dt {{AlAl) - {Al){Al)) = 0{e^) G Bn) (77) 

if (^^i^ja^ja) ^ (^ji ) (^l2 ) (^?3 ) (with the same accuracy) at t = 0. Similarly, in terms of PDF's 

9^ (P^^ (^^^ ' ) - P^f (s.i )P^f is,,)) = 0{e') in , h G Siv) (78) 

if^ff«..J^n.'*-''*«'«-)-^:'('*-)^'(--)^t'(««)^'(--)^*-0- Herep(;f,3,^(.,„.,„.,3,s.J,jff(s.. 
and Sj) are the four-mode, two-mode and one-mode PDF's obtained from P by integrating out all but 3,2 or 1 
arguments respectively. One can see that, with a accuracy, the Fourier modes will remain independent of each other 
in any pair over the nonlinear time if they were independent in every triplet at i = 0. 

Similarly, one can show that the modes will remain independent over the nonlinear time in any subset of Af < A*' 
modes with accuracy M/N (and e^) if they were initially independent in every subset of size M + 1. Namely 

(ji,j2,...,jM G Biv) (79) 

:f p(M+l,a) _ p{a) p(a) p{a) „j- j. _ n 

" ^Jl,i2,...,iAf + l - ^jl ^32 ■ ■ ■ ^JM + l at t - U. 

Mismatch 0{M/N) arises from some terms in the ZS equation with coinciding indices j. For M = 2 there is only 
one such term in the A-sum and, therefore, the corresponding error is 0(1/A) which is much less than O(e^) (due 
to the order of the limits in TV and e). However, the number of such terms grows as M and the error accumulates to 
0{M/N) which can greatly exceed O(e^) for sufficiently large M. 

We see that the accuracy with which the modes remain independent in a subset is worse for larger subsets and that 
the independence property is completely lost for subsets approaching in size the entire set, M ^ N . One should not 
worry too much about this loss because N is the biggest parameter in the problem (size of the box) and the modes 
will be independent in all M-subsets no matter how large. Thus, the statistical objects involving any finite number 
of modes are factorisable as products of the one-mode objects and, therefore, the WT theory reduces to considering 
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the one-mode objects. This results explains why we re-defined RPA in its relaxed "essential RPA" form. Indeed, in 
this form RPA is sufficient for the WT closure and, on the other hand, it remains valid over the nonlinear time. In 
particular, only property 1771 is needed, as far as the amplitude statistics is concerned, for deriving the 3- wave kinetic 
equation, and this fact validates this equation and all of its solutions, including the KZ spectrum which plays an 
important role in WT. 

The situation were modes can be considered as independent when taken in relatively small sets but should be treated 
as dependent in the context of much larger sets is not so unusual in physics. Consider for example a distribution of 
electrons and ions in plasma. The full A*'-particle distribution function in this case satisfies the Liouville equation which 
is, in general, not a separable equation. In other words, the A'^-particIe distribution function cannot be written as a 
product of A*' one-particle distribution functions. However, an M-particle distribution can indeed be represented as a 
product of Af one-particle distributions if M <^ No where Nd is the number of particles in the Debye sphere. We see 
an interesting transition from a an individual to collective behaviour when the number of particles approaches Nd ■ In 
the special case of the one-particle function we have here the famous mean-field Vlasov equation which is valid up to 
0{1/Nd) corrections (representing particle collisions). 

11 One- mode statistics 

We have established above that the one-point statistics is at the heart of the WT theory. All one-point statistical 
objects can be derived from the one-point amplitude generating function. 



which can be obtained from the A^-point Z by taking all /i's and all A's, except for \j , equal to zero. Substituting such 
values to (16411 (for the three- wave case) or to 1721 1 (four- wave) we get the following equation for Za, 




(80) 



where for the three-wave case we have: 




(82) 



(81) 



and in the four-wave case: 




(84) 



(83) 



Here we introduced the wave-action spectrum. 




{Afy. 



(85) 



(86) 



which, for p = 1 gives the standard kinetic equation. 



(87) 



First-order PDE 181)11 can be easily solved by the method of characteristics. Its steady state solution is 



Z = 



1 



(88) 



1 — Arifc 
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which corresponds to the Gaussian values of momenta 

M<'''=p!<. (89) 

However, these solutions are invalid at small A and high p's because large amplitudes s = \a\^ , for which nonlinearity is 
not weak, strongly contribute in these cases. Because of the integral nature of definitions of and Z with respect 

to the s = |a|^, the ranges of amplitudes where WT is applicable are mixed with, and contaminated by, the regions 
where WT fails. Thus, to clearly separate these regions it is better to work with quantities which are local in s = |ap, 
in particular the one-mode probability distribution Pa. Equation for the one-mode PDF can be obtained by applying 
the inverse Laplace transform to 18011 . This gives; 



dPg ^ dF 
dt dsi 



0, (90) 



with F is a probability flux in the s-space. 



SP 

F^-sAiPa + m^)- (91) 

OS J 



Let us consider the steady state solutions. Pa = so that 

-s(7P + rjdsP) = F = const. (92) 
Note that in the steady state "//rj = n which follows from kinetic equation I87II . The general solution to 1921 1 is 

P ~ Phom "4" Ppart (93) 

where 

Phom = const exp {—s/n) (94) 
is the general solution to the homogeneous equation (corresponding to F = 0) and Ppart is a particular solution, 

Ppart = -iF/ri)Ei{s/n) exp (-s/n), (95) 

where Ei[x) is the integral exponential function. 

At the tail of the PDF, s ^ n^, the solution can be represented as series in 1/s, 



F nF 

Ppart = -r-r^ + ■■■■ (96) 



Thus, the leading order asymptotics of the finite-flux solution is 1/s which decays much slower than the exponential 
(Rayleigh) part Phom and, therefore, describes strong intermittency. 

Note that if the weakly nonlinearity assumption was valid uniformly to s = cx3 then we had to put _F = to ensure 
positivity of P and the convergence of its normalisation, J P ds = 1. In this case P — Phom = n exp (—s/n) which 
is a pure Rayleigh distribution corresponding to the Gaussian wave fleld. However, WT approach fails for for the 
amplitudes s > Sni for which the nonlinear time is of the same order or less than the linear wave period and, therefore, 
we can expect a cut-off of P{s) at s = s„i. Estimate for the value of s„i can be obtained from the dynamical equation 
11411 by balancing the linear and nonlinear terms and assuming that if the wave amplitude at some k happened to be 
of the critical value s„i then it will also be of similar value for a range of fc's of width k (i.e. the fe-modes are strongly 
correlated when the amplitude is close to critical). This gives: 

s„, = Lo/eWk'' (97) 

This cutoff can be viewed as a wavebreaking process which does not allow wave amplitudes to exceed their critical 
value, P(s) = for s > s„i . Now the normalisation condition can be satisfied for the flnite-flux solutions. Note that in 
order for our analysis to give the correct description of the PDF tail, the nonlinearity must remain weak for the tail, 
which means that the breakdown happens far from the PDF core s <C s„; . When s„; ~ n one has a strong breakdown 
predicted in |32| which is hard to describe rigorously due to strong nonlinearity. 
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Depending on the position in the wavenumber space, the flux F can be either positive or negative. As we discussed 
above, -F < results in an enhanced probabihty of large wave amphtudes with respect to the Gaussian fields. Positive 
F mean depleted probability and correspond to the wavebreaking value s„i which is closer to the PDF core. When s„i 
gets into the core, s„; ~ n one reaches the wavenumbers at which the breakdown is strong, i.e. of the kind considered 
in |32| . Consider for example the water surface gravity waves. Analysis of 1321 predicts strong breakdown in the 
high-fc part of the energy cascade range. According to our picture, these high wavenumbers correspond to the highest 
positive values of F and, therefore, the most depleted PDF tails with respect to the Gaussian distribution. When one 
moves away from this region toward lower fc's, the value of F gets smaller and, eventually, changes the sign leading 
to enhanced PDF tails at low fc's. This picture is confirmed by the direct numerical simulations of the water surface 
equations reported in |23| the results of which are shown in figure Q 

This figure shows that at a high k the PDF tail is depleted with respect to the Rayleigh distribution, whereas at 
a lower k it is enhanced which corresponds to intermittency at this scale. Similar conclusion that the gravity wave 
turbulence is intermittent at low rather than high wavenumbers was reached on the basis of numerical simulations in 
|33| . To understand the flux reversal leading to intermittency appears in the one-mode statistics, one has to consider 
fluxes in the multi-mode phase space which will be done in the next section. 

12 Intermittency and the multi-mode probability vortex. 

In the previous section, we established that one-mode PDF's can deviate from the Rayleigh distributions if the flux 
of probability in the amplitude space is not equal to zero. However, in the full A'^-mode amplitude space, the flux 
lines cannot originate or terminate, i.e. there the probability "sources" and "sinks" are impossible, see 16911 . Even 
adding forcing or dissipation into the dynamical equations does not change this fact because this can only modify the 
expression for the flux (see the Appendix) but it cannot change the PDF continuity equation 1651 . Thus, presence 
of the finite flux for the one-point PDF's corresponds to deviation of the flux lines from the straight lines in the A'^- 
mode amplitude space. The global structure of such a solution in the A-mode space corresponds to a A dimensional 
probability vortex. This probability vortex is illustrated in figure ? which sketches its projection onto a a 2D plane 
corresponding to one low-wavenumber and one high-wavenumber amplitudes. Taking ID sections of this vortex one 
observes a positive one-mode flux at high k and a negative one-mode flux at low k, in accordance with the numerical 
observations of flgure|5| We should say, however, that existence of the probability vortex solutions, although consistent 
with numerics, remains hypothetical and further work needs to be done to find solutions of 16911 with non-zero curl. 

13 Discussion. 

In this paper, we reviewed recent work in the field of Wave Turbulence devoted to study non-Gaussian aspects of 
the wave statistics, intermittency, validation of the phase and amplitude randomness, higher spectral moments and 
fiuctuations. We also presented some new results, particularly derivation of the analog of the Peierls-Brout-Prigogine 
equation for the four-wave systems. The wavefileds we dealt with are, generally, characterised by non-decaying correla- 
tions along certain directions in the coordinate space. These fields are typical for WT because, due to weak nonlinearity, 
wavepackets preserve identity over long distances. One of the most common examples of such long-correlated fields 
is given by the typical initial condition in numerical simulations where the phases are random but the amplitudes 
are chosen to be deterministic. We showed that wavefields can develop enhanced probabilities of high amplitudes at 
some wavenumbers which corresponds to intermittency. Simultaneously, at other wavenumbers, the probability of high 
amplitudes can be depleted with respect to Gaussian statistics. We showed that both PDF tail enhancement and its 
depletion related to presence of a probability flux in the amplitude space (which is positive for depletion and negative 
for the enhancement). We speculated that the A^- dimensional space of A' amplitudes, these fluxes correspond to an 
A-dimensional probability vortex. We argued that presence of such vortex is prompted by non-existence of a zero- 
amplitude-flux solution corresponding to the KZ spectrum with de-correlated amplitudes. Finding such a probability 
vortex solution analytically remains a task for future. 
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14 Appendix: Wave Turbulence with sources and sinks 



One of the central discoveries in Wave Turbulence was the power-law Kolmogorov-Zakharov (KZ) spectrum, n''^ ~ fc", 
which realise themselves in presence of the energy sources and sinks separated by a large inertial range of scales. The 
exponent v depends on the scaling properties of the interaction coefficient and the frequency. Most of the previous 
WT literature is devoted to study of KZ spectra and a good review of these works can be found in We are not 
going review these studies here, but instead we are going to find out how adding such energy sources and sinks will 
this modify the evolution equations for the statistics. Instead of the Hamiltonian equation Q let us consider 

O'hC 

ici = ■^+*7iC!, (98) 

where 7; describes sources and sinks of the energy, e.g. due to instability and viscosity respectively. Easy to see that 
this linear term will not change the structure of the A'^-mode PDF equation 1651 but it will lead to re-definition of the 
fiux: 

J'.-^J'j-jjS.P. (99) 
In the one-mode equations, this simply means renormalisation 

7^7-7. (100) 

Thus we arrive at a simple message that the energy sources and sinks do not produce any "sources" or "sinks" for the 
fiux of probability. 

In the inertial range, there is no fiux modification and one can easily find F — solution of 19011 for the one-mode 
PDF, 

P,('^) = (l/nf)exp(-s,/nf), (101) 

where n*^ is the KZ spectrum (solving the kinetic equation in the inertial range). However, it is easy to check by 
substitution that the product of such one-mode PDF's, V = Tij{l/nj^) exp{—Sj /uj^), is not an exact solution to the 
multi-mode equation 1651 . Thus, there have to be corrections to this expression related either to a finite flux, an 
amplitude correlation or both. 
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Figure 1: One-mode amplitude PDF's at two wavenumbers: upper curve at ki and the lower curve at and k2 such that 
fci > fc2. Dashed line corresponds to the Rayleigh distribution. 




Figure 2: Projection of the probability flux vortex on a (si, S2) plane where si and S2 are the amplitudes at wavenumbers 
fci and k2 such that fci > /c2- 
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